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^ ■ ABSTRACT 

o 

I/"") ' We solve the Boltzmann equation for cosmological neutrinos around the epoch of 

the electron-positron annihilation in order to verify the freeze-out approximation and 
to compute accurately the cosmological neutrino distribution function. We find the 
radiation energy density to be about 0.3% higher than predicted by the freeze-out 

ON . approximation. As a result, the spectrum of the Cosmic Microwave Background 

anisotropics changes by ~ 0.3 — 0.5%, depending on the angular scale, and the 
amplitude of the mass fluctuations on scales below about 100 h~ l Mpc decreases by 
about 0.2-0.3%. 
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Q ■ 1. Introduction 

Anisotropics in the Cosmic Microwave Background (CMB) provide a powerful tool to probe 
the cosmological parameters (Hu, Sugiyama & Silk 1997). The results of four-year work of the 
COBE satellite (Bennett et al. 1996) allow us to determine the power spectrum of the CMB 
anisotropies to an accuracy of 7%, but with relatively poor angular resolution, #fwhm ~ 7°. New 
planned satellite missions, MAP (Bennett et al. 1997) and Planck (Bersanelli et al. 1997), will 
achieve an accuracy of better than 1% in the power spectrum with the sub-degree resolution 
(Bond, Efstathiou & Tegmark 1997; Zaldarriaga, Spergel & Seljak 1997). In turn, we need to 
bring the theoretical models to the same level of accuracy. 

In the standard cosmological model, after the epoch of the Big Bang Nucleosynthesis the 
relativistic particles include photons and the three species of neutrinos (Kolb & Turner 1990; 
Peebles 1993). While the abundance of photons is directly measured from the CMB observations, 
the abundance of primordial neutrinos can only be assessed theoretically. The standard way to 
perform such a calculation is to use the so called freeze-out approximation, which assumes that 
neutrinos decouple instantaneously from the rest of the universe at temperature of about 4 MeV. 
Then the distribution function of all three neutrino species retains the Fermi-Dirac form with the 
only parameter, the neutrino temperature, uniquely tied to the observed CMB temperature. 

However, even after decoupling the high-energy neutrinos still interact, albeit slowly, with the 
electron-positron plasma, contrary to the basic assumption of the freeze-out approximation. This 



interaction leads to some of the photon energy being transferred into neutrinos. But because it is 
the photon energy density that is directly measured, the total energy density of the universe in the 
relativistic species relative to the energy density in photons (which is measured observationally 
to about 0.3% accuracy) will be somewhat higher than the one predicted by the freeze-out 
approximation . 

Several previous attempts have been made to compute cosmological neutrino decoupling in 
greater detail, though still assuming that neutrino distribution functions have Maxwellian form 
(Dicus et al. 1982; Herrera & Hacyan 1989; Raha & Mitra 1991; Dolgov & Fukugita 1992), with the 
most comprehensive study given by Dodelson & Turner (1992). Recently, two more papers have 
addressed this problem with the full account for the Fermi-Dirac form of the neutrino distribution 
functions (Hannestad & Madsen 1995; Dolgov, Hansen, & Semikoz 1997). Both of these studies, 
however, have not achieved the desired level of accuracy of the numerical calculation (about 10 , 
which is equivalent to a 1% accuracy in a 1% correction to the freeze-out approximation). The 
problem is complex: solution of the full Boltzmann equation in three dimensions is at the very 
edge of modern computing capabilities. As a result, the previous calculations have only been 
able to cover slightly more that two decades in the neutrino momentum, which is insufficient to 
compute an asymptotic behavior of the neutrino distribution function. 

In this paper we complete the calculation of cosmological neutrino decoupling using an 
extensive calculation on a parallel supercomputer, placing special emphasis on achieving complete 
numerical convergence, and covering over seven decades in the neutrino momentum. 

The paper is composed in the following way. We derive and solve the Boltzmann equations 
for all three neutrino species in §y. In §||, we briefly touch upon the relevant numerical issues, 
relegating the details of our numerical method to Appendix. Finally in §||, we present our results 
for the neutrino energy density and compare them to the freeze-out approximation. We also 
obtain accurate distribution functions for the cosmological neutrinos. 



2. Neutrino Kinetics in the Expanding Universe 

The Boltzmann equation for neutrinos in the expanding universe is (Kolb & Turner 1990): 

E -f - H « 2 ik =c[a (1) 

where f u (q,t) is the neutrino distribution function, E v is the neutrino energy (E u = q since the 
neutrino mass, even if it exists, is assumed to be much smaller than our characteristic energy 
scale, ~ MeV), and C[f u ] is the collisional integral. Hereafter we use units in which h = c = 1. In 
the case of neutrinos interacting with the electron-positron pairs and other neutrino species via 
annihilation and scattering reactions, the collisional integral is (Hannestad & Madsen 1995) 

1 f d 3 P2 d 3 p3 <i 3 p 4 
2(2^f J '2E^'2El'2El 

where V\ = Q, fi = f v , 
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A(/i, h, h,h) = AA(1 - h)(l - h) - /iAU - A)(i - A), 



M 2 is the matrix element squared and summed over initial and final spin states, pi are 
four-momenta of the incoming (1,2) and outgoing (3,4) particles, and the sum is taken over all of 
the reactions involving f\. 

The list of all neutrino reactions is presented in Table [l], along with the respective matrix 
elements (Hannestad & Madsen 1995). Indecies i, j, k run over electron, muon, and tau neutrino, 
with the exception that j ^ i. The factor Gf is the Fermi coupling constant, and coefficients 
Cy and C A for different types of neutrinos are given by the following equations (for example, 
Kaminker et al. 1992): 

C v {y e ) = 2sin 2 e w + -, C A (v e ) = -, 
Cviy^Vr) = 2sin 2 9vi/--, C A (u^,u T ) = --, 

where @w is the Weinberg angle, and we adopt sin 2 Q\y = 0.23. 
Quantities Qi are defined as follows: 
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where m is the electron mass. As has been shown by Hannestad & Madsen (1995), integrals over 
<i 3 £>4 and over angles in d?p2 and d s p 3 can be computed analytically, yielding 

C V»\ = E 2(^)5 / 1^ It' A(/l ' /2 ' /3 ' W&'to'ri- (4) 

In order to minimize the possibility of an error in the complicated factors F(pi,p2,p 3 ), we have 
used Mathematica software package to perform the calculations. The resultant expressions are too 
large to be presented here, but the original Mathematica script and the FORTRAN source code 
are available upon request. 

Table 1: Neutrino Reactions 

Reaction M 2 

Vi + Pi ^ e~ + e+ 32G 2 F [(C v + C A ) 2 Q 3 + (C v - C A ) 2 Q 2 + (C 2 V - C\)Q A ] 

Vi + e" - vi + e" 2>2G\ [{C v + C A ) 2 Q l + (C v - C A ) 2 Q 3 - {C 2 V - C 2 A )Q 5 ] 

i« + e+ - i/< + e+ 32G 2 , [{C v + C A ) 2 Q 3 + (C v - C A ) 2 Q! - (C v - C%)Q 5 ] 

Vi + Vi^ Vi + Vi 128G 2 F Q 3 

Vi + v-i^ Vj + i?j 32G F Q 3 

Vi + Vj —> Vi + i?j 32G 2 F Q 3 

fi + fk^ n + v k 32G 2 F Q\ 



In order to calculate the evolution of the neutrino distribution functions, /„(C)> we need to 

include the equations describing the evolution of the scale factor and the energy density of the 

universe: 

da __ /8ttG \ l l 2 

~dt ~ \~ 

and 
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where p and p are the energy density and the pressure, respectively: 
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where A = m/T and the functions C n {\) are defined as 

C n {\) 



Vx 2 — X 2 x 2n+1 dx 
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Here p u is the energy density in the three neutrino species, 
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We also note that since the coefficients Cy and Ca are the same for muon and tau neutrino, their 
distribution functions are equal. 



3. Numerical Issues 

Equation (||), along with equations (a) and (|6|), can now be integrated numerically for each of 
the neutrino species. In order to eliminate the derivative with respect to the neutrino momentum 
in equation (|l|), we employ the comoving momentum Q = qa. We lay out the neutrino distribution 
function on a logarithmically spaced mesh in the range 10 -5 ' 5 < q/T < 10 with 40 points 
per decade (289 points altogether). By appropriately changing the limits and sampling of the 
momentum mesh, we have verified that such a discretization offers a fully convergent solution to 
an accuracy of better than 10 -4 . After this procedure, equation (|I|) becomes a system of coupled 
ordinary differential equations. We begin the integration at T = lOMeV and carry it out to 
T = 10 -3 MeV, at which point the desired precision is achieved. The relative accuracy of the 
integration at each time step is set to 10~ 7 . 

The resultant ordinary differential equations are stiff and present a considerable computational 
challenge. Standard methods require computing the full Jacobian, which is virtually impossible 
for our fairly complicated system of equations. For the purpose of this calculation, we develop 
a special numerical scheme, presented in the Appendix, which can handle stiff equations of the 
Bolztmann type and does not require computing the full Jacobian. Our scheme is more efficient 
by a factor of 20 to 60 than the standard fifth order adaptive Runge-Kutta method. 



Table 2: Main Results 



Quantity Exact value Freeze-out approximation 
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Finally, we note that the most time consuming part, calculating the collisional integral, can 
be done very efficiently in parallel. Our final computation has been performed on the NCSA Power 
Challenge Array with 12 R10000 processors and has consumed about 200 processor-hours. 



4. Results and Discussion 

The main outcome of our calculations is the number density and the energy density of 
all three neutrino species at the current epoch relative to those of photons. The results are 
presented in Table ^. For comparison, we also give the respective numbers computed in the 
freeze-out approximation. All values are accurate to the last decimal place shown. The most 
important quantity, the total radiation energy density in the universe, differs from the freeze-out 
approximation by only 0.3%. 

Another way of presenting this difference is the effective number of neutrino species, N e g. We 
can rewrite the expression for the energy density of the universe as 
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where in the freeze-out approximation N e Q = 3. From Table g we obtain: 

N eS = 3.022. 

This number does not, of course, mean that there are more than 3 species of neutrinos; it is simply 
a number that should be used in the freeze-out approximation to reproduce the exact result. Since 
most of the previously obtained results and existing numerical codes are based on the freeze-out 
approximation, it is convenient to use iV e ff : one simply has to use 3.022 instead of 3.0 in every 
place in the code where the neutrino energy density is computed. 

We note here that we find a somewhat larger effect than both Hannestad & Madsen (1995) 
and Dolgov et al. (1997), who found iV e ff = 3.017 and iV e ff from 3.013 to 3.019 (depending on 
the method of calculation) respectively. We attribute this difference to the higher accuracy of 
our calculation and the lack of numerical convergence in the previous work. In particular, when 
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Fig. 1. — The fractional difference between the effective neutrino temperature and its value in the freeze- 
out approximation, as a function of the neutrino momentum, q. Solid line is for the electron neutrino 
temperature, and dashed line is for the muon and tau neutrino temperatures. 

we adopt a momentum range from q/T = 10 _1 to 10 1 ' 3 , as in Dolgov et al. (1997), we obtain 
-Wcff = 3.019, in agreement with the authors. If we further reduce the momentum range from 
q/T = 1CP ' 3 to 10 11 , as in Hannestad & Madsen (1995), we recover their result, N e g = 3.017. 

We can also characterize the final neutrino distribution function. Let us introduce the effective 
neutrino temperature, T e g, as 

Since the neutrino distribution function is not of the Fermi-Dirac form any more, T e g is a function 
of the neutrino momentum q. Figure 1 shows the deviation of the effective temperature at the 
current epoch, To ie ff) from the freeze-out approximation value, Tq :1/ = (4/ll) 1 ' 3 To i7 , as a function 
of q/T (this ratio is independent of time after electron-positron annihilation). For the high 
neutrino momenta, the effective neutrino temperature asymptotically approaches the photon 
temperature, because the high-energy neutrinos can efficiently interact with the electrons via pair 
creation even after annihilation, 

T c g -^>T for q — ► oo. 



For the low momenta, q <C T, neutrino interactions with electrons and positrons are suppressed 
by a factor q 2 and do not affect the evolution of the distribution function. However, the rate of 
the neutrino-neutrino interactions is proportional to the first power of the momentum q. Thus, in 
the limit of small momenta, 

gjUg) _ g dT eS 

at -2*, dt Kq > 

and we find that 

T e g — ► const for q — > 0. 

The value of the constant is 0.99950 2o „ for the electron neutrinos and 0.99918 Tq m for the muon 
and tau neutrinos. 

Our results have several immediate cosmological implications. First, the change of the 
radiation energy density of the universe will affect the spectrum of CMB anisotropics at about 
0.3% level just behind the first acoustic peak, I ~ 300, and at about 0.5% level at the damping 
scale, I ~ 1000 (Hu et al. 1995). If more than 0.5% accuracy is required in calculating the CMB 
anisotropies, equation (fj) should be used. 

Second, the change of the radiation energy density of the universe affects the evolution of 
linear density fluctuations on galactic and subgalactic scales. We have computed the matter 
transfer function using COSMICS package (Bertschinger 1995) and we find that for a cosmological 
model with Qq = 1, h = 0.5, and fi& = 0.05, the rms density fluctuation at 8/i _1 Mpc scale, eg, 
decreases by about 0.2%, and the rms density fluctuation on 100/i _1 kpc scale decreases by about 
0.3%. These latter changes, however, are too small to be of any interest in the foreseeable future. 

In addition, neutrino decoupling affects primordial helium production. However, Dodelson 
&; Turner (1992) showed that the net change in the primordial helium abundance is virtually 
unobservable due to cancellation of two competing effects: one is the higher expansion rate, which 
leads to a higher helium abundance, and the other is the faster neutron decay rate, which leads to 
a lower helium abundance. A small change in the neutrino number density will produce a small 
change in the massive neutrino mass density (Hannestad & Madsen 1995), but this change is again 
too small to be of any practical interest. 

Thus, we conclude that only in computing the CMB anisotropies should one need to worry 
about the accurate calculation of neutrino decoupling; otherwise, the effect is negligibly small. 
Overall, our calculations confirm the validity of the freeze-out approximation. If the accuracy of a 
few percent is sufficient, one can safely use the freeze-out approximation to compute any property 
of cosmological neutrinos. 

We are grateful to D. Yakovlev, J. Ostriker, D. Spergel, J. Madsen, B. Fields, S. Hannestad, 
and A. Dolgov for valuable comments. We thank the anonymous referee for pointing us to an 
error in the original manuscript. N. G. was supported by the UC Berkeley grant 1-443839-07427. 
Calculations were performed on the NCSA Power Challenge Array under the grant AST-960015N. 



A. Numerical Method 

In this paper we are dealing with a particular kind of an ordinary differential equation that 
can be presented in the following form: 

^ = f(y) = w(y) - k(y)y, (Al) 

where both w and k are slow functions of y, but not of t, and k is positive. This equation is stiff, 
and a numerical method which does not handle stiff equations requires a time step At such that 

kAt <C 1. 

Numerical methods that can deal with stiff equations usually have a much less stringent restriction 
on the time step, 

Ai< 1, 



<C k. 



dy 
and, because w is a slow function of y, we assume that 

dw 
dy 

However, standard techniques for stiff equations require computing the full Jacobian, 

_ dw dk 
dy dy 
In our case this quantity is very difficult to compute, because w is an integral over y and numerical 
evaluation of the integral involves nontrivial interpolation. 

We therefore proceed differently and design a numerical scheme which involves only a partial 
Jacobian, 

J=-k, 



which can be computed simultaneously with the r.h.s. of equation (Al) at no extra cost. 



Additional advantage of using the partial Jacobian J instead of the full Jacobian J is that 
for a system of equations, the partial Jacobian is a diagonal matrix, which can be inverted much 
faster than the full Jacobian, which is usually a general matrix. 

However, we cannot simply take a standard numerical scheme and replace the full Jacobian 
with the partial one, because the different orders of the numerical error will not cancel out in this 
case. Thus, we need to design a special scheme which will assure the proper cancellation of the 
numerical error up to a given order. 

The numerical scheme to update y from y = yo to y = y\ in a time interval h is constructed 
as follows: 

hf(yo) 
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1 + 7/ifc' 

hf(yo + a2igi) + c 2 igi 

1 + jhk 
hf(yo + 03101 + 03252) + C3101 + C3252 
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2/1 = Vo + hgi + b 2 g2 + h93, (A2) 



Table 3: 



Table 4: 



Constant 



Constant 


Value 


7 


0.788675134594812882251 


Q21 


1 


Q31 


0.56698729810778067662 


032 


1/4 


C21 


-1.26794919243112270647 



Third order scheme 



Second order scheme 



cai -3/2 

c 32 -1.1830127018922193234 

6i 1.3779915320718537844 

& 2 0.9553418012614795489 

b 3 2/3 



-3.183012701892219323 
-3.049038105676658006 
1.566987298107780677 
1.038675134594812883 
0.21132486540518711775 



where 7, Oj, 6j, and q are constants. The values for 7, 021, 031, 032, and C21 are given in Table pJL 

By varying the remaining constants, we construct two numerical schemes: the third order and 
the second order, respectively (Table Q). Thus, the difference between the values of y\ computed 
with the two schemes can serve as an estimate of numerical errors. Again, Mathematica was 
used to compute the values of the constants that give the cancellation of numerical errors to the 
required order. 
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